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Nonlinear dynamics of the magnetic vortex state in a circular nanodisk was studied under a perpendicular 
alternating magnetic field that excites the radial modes of the magnetic resonance. Here, we show that as the 
oscillating frequency is swept down from a frequency higher than the eigenfrequency, the amplitude of the 
radial mode is almost doubled to the amplitude at the fixed resonance frequency. This amplitude has a 
hysteresis vs. frequency sweeping direction. Our result showed that this phenomenon was due to a 
Duffing-type nonlinear resonance. Consequently, the amplitude enhancement reduced the vortex 
core-switching magnetic field to well below 10 mT. A theoretical model corresponding to the Duffing 
oscillator was developed from the Landau-Lifshitz-Gilbert equation to explore the physical origin of the 
simulation result. This work provides a new pathway for the switching of the magnetic vortex core polarity 
in future magnetic storage devices. 

A magnetic vortex state^'^ is a ground state of a magnetic nanostructure that consists of a perpendicularly 
magnetized core and in-plane curling magnetizations around the core (Fig. la). Because of its importance 
in fundamental physics, research on the vortex state is an important emerging topic in magnetism 
studies^"^ and it has a high potential for application in high-density data storage devices^"^^ A magnetic vortex 
state is energetically fourfold degenerated, which is determined by its polarity and chirality, where the polarity 
refers to the perpendicular direction of the core magnetization, pcore (with a value of + 1 for the upward direction, 
or -1 for the downward direction) and the chirality, c, refers to the curling direction of the in-plane magnetization 
(a value of + 1 or -1 for the counterclockwise or clockwise directions, respectively). Obviously, the success of a 
magnetic vortex device will critically depend on the question of how to control the vortex polarity and chirality 
effectively. Much effort has been invested recently in developing various methods for reversing the vortex polarity 
and chirality with a low magnetic field. While the chirality can be reversed easily with a weak field of ~ 50 mT (ref. 
13), the magnetic field required to reverse the vortex core is on the order of 500 mT, which is too large for practical 
use in device applications^^"^^. 

To reduce the vortex core-reversal field, an alternative approach used a dynamic field^"^°'^^'^^"^^. The chirality 
reversal field can be almost halved when a short-pulsed field is applied^^ A promising result is also reported for an 
AC oscillating magnetic field set at the vortex resonance frequency, so that the vortex excitation could assist its 
polarity reversal. A representative example of such an approach is the vortex gyration excitation, in which the 
vortex core exhibits a spiral motion as an AC magnetic field is tuned on at the gyration eigenfrequency^"^°'^^'^\ 
Core switching occurs subsequently through vortex-antivortex creation and annihilation^ as the core's moving 
speed exceeds a critical value^^. The core reversal field can be reduced in such a manner to values far below 10 mT 
(refs 6,7,19). 

However, this method contains a fundamental problem for applications. After the core reversal and turning off 
the field, the core gyration exponentially decays to its initial position. The decay radius is comparable to the lateral 
size of the sample and the relaxation takes a few hundred nanoseconds^^ This is a severe obstacle to reading the 
polarity^^. 

Recently, Wang and Dong^^ and Yoo et al}^ found a new method of vortex core flip from numerical simulation. 
They demonstrated that the vortex core polarity could be switched in a radial excitation mode^^'^^'^^'^^ by a 
perpendicular AC magnetic field. In contrast to the gyration mode-assisted switching, which involves the vortex 
core motion, the radial mode-assisted core switching involves only axial symmetric oscillations, thus preserving 
the vortex core position. Obviously, the radial mode-assisted core switching has a completely different mech- 
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Figure 1 | Magnetic vortex and core reversal induced by sweeping 
frequency, (a), Disk structure used for simulations and the initial vortex 
state. The in-plane curling magnetization directions are shown as arrows. 
At the disk center, a perpendicularly magnetized core exists, (b), The 
magnetization oscillation induced by the external magnetic field 
applied to the thickness direction (z). The field strength was 9 mT with a 
sinusoidal form; the frequency was varied with 0.2 GHz/ns rate; <mz> is 
the spatially averaged magnetization component in the z direction 
normalized by the saturation magnetization, (c), I^, the amplitude of 
<mz> oscillation, as a function of external field frequency / The 
amplitude obtained from downward sweeping (red line) can achieve 
sufficient amplitude for core polarity (pcore) reversal but from upward 
sweeping (blue line) it cannot. Black circles denote saturated amplitude per 
each frequency. This fixed frequency result also does not obtain the 
maximum amplitude for pcore switching. Snap shots show the 
magnetization states. The overlapped images represent one cycle of the 
oscillation and single images exhibit the reversal process of the core. 

anism from the gyration mode-assisted core switching. The under- 
lying mechanism of the radial mode- assisted core switching was not 
clearly shown by the simulation. The critical field obtained by the 
radial mode in these studies is of the order of 20 mT (ref 17), larger 
than the gyration mode- assisted core reversal. 

In this work, we studied the underlying mechanism of the radial 
mode oscillation and outlined a new pathway to reduce the core- 
switching field further down to the mT range, which was more com- 
parable to the critical field of the gyration -assisted core switching. In 
addition to micromagnetic simulations^^, we also established a 
dynamical equation for the radial mode oscillation from the 
Landau-Lifshitz-Gilbert (LLG) equation^^. This equation clearly 
explores the nonlinear behavior of the radial mode and the critical 
field reduction. 

For direct comparison of the critical field reduction, the simu- 
lation structure was set as described by Yoo et aV^ (Fig. la). 
According to previous studies, the radial modes are classified by 
the node number n (refs 17,24). The first mode has one node, the 
vortex core, which means that the magnetization does not oscillate 
temporally at the vortex core, but the other parts almost uniformly 
oscillate. The second mode has two nodes; one is the vortex core and 



the other a concentric circle. Yoo et aU^ studied the resonance fre- 
quency of the individual radial mode and obtained the eigenfrequen- 
cies with the same sample structure as in this study: 10.7 GHz for the 
first mode (n = 1), 15.2 GHz for the second mode (n = 2), and 
20.7 GHz for the third mode {n =3). They also showed the vortex 
core polarity reversal using the first mode with an oscillating external 
field of 20 mT. 

To reduce the radial mode-induced critical field below 10 mT, we 
stimulated the first mode of the radial oscillation with a different 
method; that is, sweeping of the external field frequency. The field 
was sinusoidal with amplitude of 9 mT and the field frequency /was 
slowly varied from 14.0 to 6.0 GHz over 40 ns. Figure lb shows the 
magnetization oscillation during frequency sweeping with time. The 
normalized magnetization along the thickness direction and the 
external magnetic field, H^, were plotted together. The term <m^> 
means the spatial average over the entire disk. The magnetization 
oscillation has the same frequency despite the phase difference. 
From this oscillation, we can get the oscillation amplitude of magnet- 
ization, 4, in the thickness direction, which is half the difference 
between the nearest maximum and minimum values of the <m^> 
oscillation. After reaching an external field frequency of 6.0 GHz, the 
frequency sweeping direction was reversed and/returned to 14.0 GHz. 

In Fig. Ic, 4 is shown as a function of/. It is interesting to note that 
an external field of 9 mT can reverse the vortex core polarity. In 
downward sweeping of the frequency, the almost uniform magnet- 
ization oscillation was observed on the disk except for the core con- 
serving its width (inset of Fig. Ic). This uniform oscillation was 
maintained before 4 reached the maximum amplitude of 0.28 when 
/was 8.7 GHz. After reaching this critical amplitude, the uniform 
oscillation collapsed and converged into the disk center that gener- 
ated a breathing motion of the core. Such a breathing generated a 
strong exchange field when the core was compressed, and then core 
polarization switching occurred^^'^^. Amplitude fluctuations near 
8.5 GHz and 10.5 GHz are transition effects discussed below. In 
contrast to downward sweeping, the upward frequency sweeping 
did not reach the amplitude of 0.28, so the vortex maintained its 
polarity. This means that one cycle of frequency sweeping generated 
one core reversal. It is notable that the amplitude obtained with the 
frxed-field frequencies was the same as the upward sweeping. The 
fixed-field frequency amplitudes were determined by amplitude sat- 
uration after turning on the external oscillating field. To reverse the 
core polarity with the upward sweeping oscillation and fixed- 
frequency oscillation, a larger field was required for achieving the 
sufficient oscillation amplitude. From this sweeping frequency simu- 
lation, it was verified that the critical field was reduced to below 
10 mT and this reduction was only observed in downward sweeping 
because of the hysteresis behavior of the frequency. 

To study this hysteresis behavior, we constructed a simplified 
model. The magnetization oscillation of the first mode can be repre- 
sented by only two variables, 0 and (p, because with the exception of 
the core region, the magnetization oscillates almost uniformly^^. 
Figure 2a shows the definition of 0 and (p in the core-free model. 
The angle 0 represents the magnetization tilting toward the normal 
direction and (p denotes tilting in the radial direction. The initial 
magnetization is then described as ^ = 0 and g? = 0. Note that the 
initial magnetization of the vortex state could have two values of cp, 0 
and 71, because of the chirality of the vortex. The two states are 
energetically equivalent, so in this paper we only considered (p = 0 
as the initial state. Using these two variables, the magnetization state 
of each component was described as = Ms cosO sin^, = Ms 
cosO cos^, and = Ms sinO. After inserting these components into 
the LLG equation^^ and solving, we could derive the following 
coupled equations: 

g) = yMs (Nz — Nr sin^ (p) sin 9 + yH sin cot H ^ 'O, (1) 
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Figure 2 | Simplified model for the first radial mode, (a), The angle 6 is 
tilting toward the perpendicular direction and (p represents the radial 
directional tilt, (b), A spiral trajectory during relaxation; before the 
relaxation, several Gaussian pulses were applied to excite the 
magnetization oscillation. (c),(d), Relations between two coupled 
variables, 6 and ^, which were obtained from (b). 

6 = — yMsNr cos 0 sin g? cos (p — a cos Og?. (2) 

Here, the over dot means the time derivative, y is the gyromagnetic 
ratio, Nz and Nr are determined by the demagnetizing energy along the 
z and r directions, respectively, H is the external field amplitude, and 
CO (= 2nf) is the angular frequency of the external field. In equations 
(1) and (2), all the parameters are known except for and N^. To 
obtain these two parameters, we reduced equations (1) and (2) assum- 
ing that Nz > Nr because of the thin structure, and this condition also 
results in ^ ^ 1, after which simplified equations are derived: 

g? = yMsNzO + yH sin ojt + dO, (3) 
^ = — - yMsNr sin 2(p — a^. (4) 



We carried out a simulation that solved equations (3) and (4). 
Gaussian-shaped external field pulses were applied sequentially to 
the nanodisk to excite the first mode and the interval between pulses 
was tuned for precise resonant amplification^ \ We turned off the 
external field pulse and observed the relaxation process (Fig. 2b). The 
angles 0 and (p converged exponentially to zero with a spiral traject- 
ory; 0 was obtained from sin^ = <mz(t)> — <mz(0)>, where, t is 
time and g) was determined by the averaged value over the entire disk. 
From this spiral relaxation, g) was plotted as a function of 0 (Fig. 2c) 
and sin2^ was plotted with 0 (Fig. 2d). The damping effect was 
negligible because a = 0.01 (^1) and during relaxation, H = 0. 
Thus, two linear relations (Fig. 2c and 2d) determined the values 
Nz = 0.85 and Nr = 0.14. These values are reasonable because Nz 
+ AT^ ~ 1 for the first mode was composed with almost uniform 
magnetization. 

Next, we solved equations (3) and (4) to obtain the oscillation 
amplitudes versus the external field frequencies. The time derivative 
of the equations and the elimination of the 0 and 0 terms yielded the 
following equation of motion of g): 

(1 + a^)(p + -f sin 2g) + C(iV, + Nr)g) 

2 (5) 

— C(2Nr sin^ g))g) — yHco cos cot = Q. 




f(GHz) f(GHz) 

Figure 3 | Comparison between solutions of the Duffing equation and 
simulations. The external field was set to be 8 mT. (a), Solutions of 
equation (6), the Duffing equation. The plane pendulum asymptotically 
corresponds to the Duffing oscillator (inset). The solution comprises a 
stable and an unstable line and where the stable line is disconnected, the 
amplitude and the phase d show drastic jumps, g)^ has the maximum 
value ^max at the zero phase = 0). (b), Micromagnetic simulation results 
with the same external field show similar hysteresis behavior. 

Here, ojq = yMs\/NzNr= 2n X 10.7 GHz and C = ayMs. If we 
neglect the damping term and the external field term, equation (5) 
exactly corresponds to the simple plane pendulum, O + cOq sin O = 0 
with O = 2g), so it was natural for the main dynamics of the radial 
mode of the vortex to be the same as that of the plane pendulum. 
Directly solving equation (5) was not easy, so we simplified equation 
(5). Firstly, we compared the third and the fourth terms of equation 
(5), both terms have g) dependence, and we neglected the fourth term 
which is much smaller than the third term. Then, we expanded sin 2g) 
in the second term by using a Taylor expansion such as sin x = x — 
x^/3\+x^/5\. . . and we adopted up to the third order of this expansion. 
These simplifications produced the following equation. 

(p-\-colg)(l-\- [^g)^)-\-dg) — F CO cos cot = 0. (6) 

This is the well-known Duffing oscillator equation^^'^^ that describes 
a nonlinear oscillation induced by the g? variation -dependent spring 
constant. If j5 = 0, equation (6) becomes a harmonic oscillator. In this 
study, = -2/3 due to the Taylor expansion of sin 2g). Other para- 
meters were determined similarly: the dissipative constant d = 
(xyMsiNz + Nr) and the external force F = yH. It is well known that 
we can solve the Duffing equation by using the van der Pol 
transformation^^'^°. 

A solution of equation (6), the frequency- dependent amplitude g)o, 
exhibits hysteresis behavior (Fig. 3a). We set H to be 8 mT. In con- 
trast to a harmonic oscillator showing a symmetric resonant peak, 
the Duffing oscillator exhibits an asymmetric resonance peak and the 
curve is divided into stable and unstable solutions that are well 
known as the foldover effect. If we apply the external field with a 
starting frequency of 14 GHz and monotonically reduce the fre- 
quency, the oscillation amplitude would follow the stable solution 
line until it meets the maximum amplitude point ^max- After passing 
the maximum point, the amplitude was drastically reduced to follow 
the stable line. A similar behavior was also expected for the case of 
increasing frequency, but the amplitude did not reach ^^lax because 
the stable line was connected up to 10 GHz. This hysteresis property 
was also shown as the phase 3 variation of g? oscillation with respect 
to the field oscillation (Fig. 3a). The micromagnetic simulation result 
showed almost the same hysteresis behavior (Fig. 3b) with the 
Duffing oscillator, but it showed fluctuation when the amplitude 
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and the phase jumped from one stable solution to the other stable 
line. Note that in Fig. 3 the amplitude is shown with (p. It differs from 
Fig. Ic, which was plotted with 4 sin^o)» but the maximum values 
of (Pq and can be converted through the relation Nz sin^ ^^ax — 
sin^ ^max because the demagnetization energy is alternatively trans- 
ferred between the z and r directions. 

We observed the Duffing- type oscillation and critical field reduc- 
tion in other radial modes. Figure 4 presents the minimum field for 
core reversal, H^, as a function of frequency/. To obtain the critical 
field with fixed frequency, we fixed the external field frequency and 
varied the amplitude of the field from 0 with a 10 mT/ns rate. 
Frequency dependent critical field had three local minima and each 
minimum corresponded to the first, second, and third mode. The 
obtained values for modes with fixed frequency were respectively 
20 mT for w = 1, 78 mT iox n = 2 mode, and 127 mT for the n = 3 
mode. The critical fields were also determined by the sweeping fre- 
quency method. For the second mode, the frequency started from 
17 GHz with —0.2 GHz/ns sweeping rate; for the third mode, it 
started from 24 GHz with the same rate. Then, the critical field 
was 8.5 mT for n = 1, 37 mT for n = 2, and 76 mT for n = 3, which 
were almost half the values obtained with fixed frequency. The 
sweeping rate can change the critical field but the variation is small. 
For example, n = \ mode. He = 8.3 mT with 0.1 GHz/ns rate and 
8.6 mT with 0.4 GHz/ns rate. The critical field is obtained with 
0.1 mT field interval for n = 1 and 1 mT interval for n = 2,3. We 
confirmed that these reduced critical fields were not achieved 
through the upward sweeping of the field frequency because of the 
Duffing-type nonlinear oscillation. In addition, further reduction of 
the critical field was achieved by using a square wave type external 
field. For n = 1 mode, we obtained 6.3 mT for the vortex core 
reversal. 

We tested the scalability of the radial mode-induced core reversal. 
When the radius of the disk was 120 nm, the critical field obtained by 
the frequency sweeping method was 9.3 mT. The core of a disk with 
radius 250 nm reverses its polarity with 12 mT external field. 
Increasing the radius, the critical field also increases. This scalability 
is an important property for developing data storage devices. 
Contrary to the radial mode-induced polarity switching, the critical 
field with the gyration-induced polarity switching exhibits inverse 
radius dependence as well as the chirality reversal^ ^. 

Finally, we point out the chaotic behavior and the phase commen- 
surability in the radial mode oscillation for further studies. Petit- 




f{GHz) 

Figure 4 | Critical fields for the core reversal with respect to the field 
frequency £ The fixed /results were obtained by amplitude sweeping 
simulations and the results show three local minima that correspond to 
each radial mode (n = 1,2,3);(#) denotes the critical fields determined 
by the sweeping frequency method; dashed lines represent the maximum 
amplitude-frequency points in the sweeping method. 



Watelot et al. observed the chaos and phase-locking phenomenon 
in the vortex gyration with the core reversaP\ We observed similar 
behavior in radial mode oscillation. It is expected that a nonlinear 
oscillator with a sufficiently large driving force would exhibit chaotic 
motion. We confirmed this chaotic behavior in the radial mode of the 
vortex. When the oscillating field strength was smaller than H^, a plot 
of the variable with respect to its time derivative, for example < > 
versus <m^>, showed a circular trajectory. But when the field was 
larger than H^, this plot becomes complex in the phase space, which 
manifests its chaotic behavior. Figure 5 shows examples of the chaos 
in the radial mode. The frequency was fixed at 13.5 GHz. When H = 
60 mT < He (Fig. 5a), it showed a closed circular trajectory, but when 
H = 90 mT > He the trajectory was not closed (Fig. 5b). Further 
increases in the field resulted in closed trajectories. However, the 
trajectories were not a simple circle. To close the trajectory, 14 cycles 
of field oscillation are needed (Fig. 5c) and during these 14 cycles, the 
core reversed four times. In the case of H = 120 mT, core reversal 
occurred twice in five field oscillations (Fig. 5d), implying that the 
core reversal rate was related to the chaotic behavior. Thus, to 
describe the radial mode of vortex including its chaotic behavior, 
the core polarity- related term^^ is needed in the motion equation. 

In summary, we studied the nonlinear resonance of the radial 
mode of the vortex and found that the oscillation mode correspond- 
ing to the Duffing-type nonlinear oscillator exhibited a hysteresis 
behavior with respect to the external field frequency. Through the 
hysteresis effect, we can achieve hidden amplitude that is almost 
double that obtained with fixed field frequency and this amplitude 
multiplication effect reduces the critical field below 10 mT. In addi- 
tion, we pointed out the chaotic behavior of the radial mode for 
further studies. We think that to complete the study on vortex 
dynamics, it is timely to start research on the nonlinear behavior 
in radial modes, as well as in other oscillations of the magnetic vortex. 

Methods 

Micromagnetic simulations. We performed a micromagnetic simulation using the 
OOMMF code^*^ for numerical calculation of the LLG equation. The simulation 
structure is a circular disk with a diameter of 160 nm and 7 nm thickness, t^, as shown 
in Fig. la. The structure is divided into 2X2X7 nm^ unit cells. Uniform 
magnetization was assumed along the thickness direction. The material parameters 
were chosen to resemble a typical permalloy with no crystal anisotropy. The 
saturation magnetization, M5, was 8.6 X 10^ A/m, the exchange stiffness. A, was 1.3 X 
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Figure 5 | Chaotic behaviors of the radial mode oscillation on the phase 
space. The external field frequency was 13.5 GHz. Phase space plots of 
<mz> are shown in (a-d). We used different colors, red and blue, to 
distinguish the core polarization, (e-h), Poincare sections with zero phase 
(3 = 0). To obtain a steady state trajectory, simulation data from 10 to 
20 ns was plotted, except for (b) and (f). We plotted (b) and (f) with the 
data from 10 to 500 ns. (a), When the applied field strength is smaller than 
Hey the trajectory resembles a closed circle and the Poincare section is a 
single point (e). (b), The trajectory fills some area of the phase space with 
the increased field. (c),(d), Further increase of the field results in closed 
trajectories, but they do not show a simple circle. The numbers of points in 
the Poincare section represent the number of field cycles needed to close 
the trajectories. 
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10"^^ J/m, and the damping constant a was 0.01. Then, the magnetization state of the 
structure had a well-known vortex state composed of a perpendicularly magnetized 
core at the center of the sample and curling magnetizations around the core. The 
oscillating external field was applied to the thickness direction, and then the field 
mainly excited the radial mode that also had axial symmetry. 
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